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ABSTRACT 

A mass model that includes galaxies in and near the Local Group and an 
external mass in the direction of the Maffei system, with the condition from 
cosmology that protogalaxies have small peculiar velocities at high redshifts, 
allows a plausible picture for the past motion of the Large Magellanic Cloud 
relative to the Milky Way. The model also fits the proper motions of M33 and 
ICIO. 



1. Introduction 

Analyzing the past orbit of the Magellanic Clouds relative to the Milky Way is an 
interesting problem in dynamics, and it may be a useful preliminary exercise for analyses 
of ongoing advances in measurements of galaxy distances and proper motions that will test 
ideas about the evolving mass distribution around galaxies. Since the positions, velocities and 
accelerations of the Magellanic Clouds are fairly well constrained at low redshifts (Kallivayalil 
et al. 2006, Piatek, Pryor & Olszewsk 2008), the challenge here is to understand how the 
Clouds were directed to their present paths by interactions among newly forming galaxies 
at high redshift. This analysis of the motion of the Large Magellanic Cloud (LMC) around 
the Milky Way (MW) focuses on the influence of galaxies in and near the Local Group. The 
Small Cloud is taken to be an unimportant perturbation at the hoped for accuracy of this 
study. 

The flrst steps to the interpretation of the measured proper motion of the Large Mag- 
ellanic Cloud considered its interaction with MW and M31 (Besla et al. 2007; Shattow & 
Loeb 2009; Kallivayahl et al. 2009). Peebles (2009, P9) added the initial condition indicated 
by the gravitational growth of structure in the standard Big Bang cosmology, that peculiar 
velocities of the protogalaxies at high redshift are growing. P9 applied this condition to a 
illustrative mass model. Here the analysis is extended to a more realistic model that takes 
account of the accepted picture of the mass structure in the Local Group (LG). Sensitivity 
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to the mass outside the LG is explored by comparing solutions with and without a large 
external mass near IC342 and Maffei 1. 

The numerical method grew out of the action approach in Peebles (1989) for dealing with 
the mixed boundary conditions that initial peculiar motions are small and present positions 
of the galaxies are consistent with what is observed. This paper introduces a method of 
solution that improves the efficiency of the numerical action method enough to suggest a 
new name, NNAM. The method of solution is described in Section [2] and the Appendix. 

Since the computation is approximate — the broadly distributed mass in a galaxy 
is represented by a single particle, and the model deals with only a few galaxies — one 
cannot expect to achieve an exact fit to the galaxy positions and motions. To deal with 
this, the approach taken here is to seek orbits that minimize a measure of fit to the 
physical parameters with assigned uncertainties that are treated as standard deviations. The 
physical parameters include the proper motions of LMC, M33 and ICIO, where the stated 
measurement uncertainties are treated as standard deviations. The circular velocity of MW 
is assigned what seems to be a reasonable measurement uncertainty. The galaxy masses 
are assigned central values and uncertainties, as in Shaya & Peel (2007), that are meant to 
represent the ranges of what are generally considered to be reasonable values. The present 
positions and redshifts of the galaxies are allowed uncertainties larger than the measurement 
errors. This may may be taken to represent an actual offset of the center of starlight from 
the effective center of the mass that is dominated by the dark matter on the outskirts. For 
example, if the dark matter 100 kpc from the center of M31 were as irregularly distributed as 
the stars (McConnachie et al. 2009) then it would be easy to imagine a substantial offset of 
the effective mass center of this galaxy from the center of its stars. But these offsets in galaxy 
positions and motions may just as well be considered a way to allow for the approximate 
nature of the model. The parameters and their assigned uncertainties for this procedure are 
discussed in Section |3l 

This approach allows considerable freedom of adjustment of parameters, but there is a 
considerable number of constraints. The results in Section |4] include discovery of just one 
form of orbits, within modest variations, that offers a reasonable fit to all the parameters, 
including the proper motions of LMC, M33 (Brunthaler et al. 2005) and ICIO (Brunthaler 
et al. 2007). I argue that the results offer a good case for the shape of the orbit of LMC 
relative to MW back to high redshift. The motion of M31 requires more work, and since 
M33 and ICIO are affected by what M31 has been doing the orbits presented here for M33 
and ICIO are best considered examples of what might have happened, to be checked by more 
ambitious models and still better galaxy distances and proper motions. 
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2. Numerical Methods 

This section reviews the dynamical model and the numerical action procedure. The 
details of the new method of solution are described in the Appendix. 

The starting idea is to reduce the evolution of the mass distribution to an N-body 
problem in which a single particle represents the mass concentrated around a galaxy or a 
tightly bound system of galaxies. This is thought to be a good approximation at low redshift 
because the galaxies in our neighborhood by and large are well separated relative to standard 
estimates of the sizes of their dark matter halos. The approximation need not be vitiated 
by merging, because a particle may trace the effective mass center and momentum of the 
two systems before the merger as well as in the merged system. Further details of this line 
of argument may be traced back through P09. 

In this analysis of the N-body problem the Friedmann-Lemaitre cosmology is taken to 
be spatially fiat with constant dark energy density. The expansion parameter a{t) satisfies 

with present value Oq = 1. Hubble's constant is H^, fl is the density parameter, and matter 
pressure is ignored. 



The action 



S 



(2) 



summed over particles i, with fixed present positions and the initial condition 



^ at a{t) 0, (3) 



gives the equation of motion 



d_yj^^j QH^Y (4) 



dt dt V 2a2 

The Cartesian coordinate label is k = 1,2,3. Physical lengths are Vi^k = ctXi^k, and gi^k = 
Qi^k/O''^ is the physical acceleration of particle i caused by the gravitational attraction of the 
other particles. Since the cosmic mean mass density is p(t) = 3QH^/{8TTGa{tY), the sec- 
ond term in parentheses on the right hand side of equation Q is the physical acceleration 
AnGpTi^k/'i- This counter term eliminates the accelerations of the particle coordinate posi- 
tions Xi^k when the coordinate positions are at rest and arranged so their mutual gravitational 
accelerations gi^k are the same as that of the background cosmology. 
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In a discrete time step approximation the equation of motion ^ is 

n_ C _ «n+l/2<^n+l/2 _ 



+ 



^n+1/2 " tn-l/2 



^ ^ 9i,j,k,n ~l~ c^^H^Xi^k;, 



(5) 



an-l/2""-l/2 ^ 
~r l^3;j,fc,n Xi,k,n—l)- 



The coordinate acceleration is here written as the sum over the pull of all other particles, and 
it will be recalled that the physical acceleration is gij^k = Qi.j.k/o^- The indices are Cartesian 
coordinates /c = 1, 2, 3, 1 < i < the particle number, and time steps 1 < n < n^; + 1. The 
expansion parameter a„+i/2 and time tn+i/2 are evaluated half-way between time steps n 
and n + 1, producing a leapfrog approximation to the equation of motion. In the action 
method the present positions Xj are given and the initial condition in equation ([s]) is 

represented by setting ai/2 = so that equation ([s]) at n = 1 is 



'^3/2'^3/2 ^3/2 / X — ^ 1 9 \ 

= Si^k,i = {xi,k,2 - H ^ fi-ij^,! + -QH^Xi^k,i ■ (6) 

CL2 — CLi ■^"^ 



When the free coordinates Xi^k,n at 1 < n < are adjusted to make the same number 
of Si^k,n all vanish then the Xi^k,n are a numerical solution in leapfrog approximation to 
the equation of motion Q for given final positions and the initial condition in 

equation One approach to Si^k,n = is to walk down the the gradient of S* to a minimum 
(Peebles 1989), but that misses maxima and saddle points. The Si^k,n are driven to zero at 
these stationary points of the action by iteratively applying the coordinate shifts 5xi^k,n from 
the solution to 

Si,k,n + ^ 6Xj^k',n' = 0. (7) 

J,k',n' ' 

Computing the Sxi^k,n by matrix inversion, as in Peebles et al. (2001), is exceedingly slow 
for an interesting number of time steps. The NNAM described in the Appendix speeds this 
up by making use of the fact that in equation ^ the matrix dSi^k,n/dxj^k',n' is nonzero only 
near the diagonal. 

This analysis of the motion of LMC also uses numerical integration of the equation of 
motion ^ forward in time from initial positions and proper peculiar velocities taken from 
a NNAM solution. These initial conditions are 

/ I \ /o ■ • ^i,k,2 ~ Xi^k,l /„N 

Xi,k.3/2 = [Xi^k,2 + 3;i,fc,lj/2, Vi,k,3/2 = 0'3/2Xi,k,3/2 = 03/2^3/2 • [o) 

0,2 ~ O'l 
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This is accurate when the time steps are uniformly spaced in the expansion parameter a{t) 
and the particle coordinates are varying in proportion to a{t), as in the growing mode of 
departure from homogeneity of an ideal fluid in linear perturbation theory. It proves to be 
a good approximation in the numerical solutions presented here, as evidenced by the close 
agreement of present positions given to NNAM and present positions and velocities from 
forward numerical integration with these initial conditions. 



3. Parameters 

3.1. Numerical parameters 

To take account of the dark matter halos of MW and LMC when LMC approaches 
the MW at low redshift their mutual interactions are modeled on a rigid truncated limiting 
isothermal sphere, with coordinate acceleration 

avl mi^Mcavl , Guimw 
5'LMC = , 5'MW = , for ax <ri = — , (9) 

where niMw and mLMC are the masses and Vc is the LMC circular velocity at the Solar 
circle. This prescription conserves momentum. It ignores the perturbation to the MW mass 
distribution by LMC, as in dynamical drag, but P9 shows that the effect of dynamical drag 
on the motion of LMC is small in the LG model used here, because LMC has approached 
at high velocity. For other particles, and for the interaction of LMC and MW at larger 
separation, the coordinate acceleration of particle i due to j is the usual inverse square law 
form, 

9id,k = Gmj ' _ ' (10) 

I I J I 

The MW mass structure is known in more detail than the model in equations ^ and ( [Io| , 
but this simplified scheme seems to be an appropriate match to the schematic nature of the 
mass model. 

The NNAM solutions have = 500 time steps, about as many as convenient compu- 
tation time allows. The numerical integration forward in time from the initial conditions in 
equation (|8| has rix = 5000 steps. In both cases the expansion parameter a„ is uniformly 
spaced. The NNAM solutions start at = 0.1, or initial redshift 

l + ^i = 10. (11) 

In the iterative relaxation to a NNAM solution to the equation of motion the coordinates 
are shifted by e5xi^k,n where 5xi^k,n is the solution to equation ([7|. At e = 1 the iteration on 
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occasion fails to converge. This is avoided by setting e = 0.1 far from a solution or where 
convergence is unusually slow. In the results presented here the orbits are relaxed to the 
equation of motion Q to near machine double precision. 

The search for plausible NNAM solutions starts from random orbits (with parameters 
uniformly distributed over the ranges of nominal standard deviations in Table 1, orbits of 
LMC, M33 and ICIO that circle present positions Xq as x = Xo + Asm{4:7rB{l — a{t))), where 
A is uniformly distributed between ±5 Mpc and B is uniformly distributed between and 1, 
and the other galaxies move linearly in a{t) to present positions from initial displacements 
uniformly distributed between ±5 Mpc in each Cartesian coordinate). If the NNAM solution 
looks promising it is used to get initial conditions from equation Q for a numerical inte- 
gration of the equation of motion forward in time, with an order of magnitude smaller time 
step. In the solutions discussed in Section |4] present positions and velocities from NNAM and 
the forward integration never differ by more than 0.2 kpc and 1 km s~^, but for security the 
latter with the smaller time steps is used. The masses, present positions, and MW circular 
velocity Vc are iteratively shifted, one at a time, to get a new solution, first from NNAM 
and then the forward integration. If a parameter shift decreases a measure of fit the shift 
is increased by 25% for the next application. If it increases the shift is halved, the sign 
reversed, the parameter reset to half the distance between original and trial values, and 
recomputed. 

The simple minimization procedure is adequate for the present purpose but I expect 
ought to be improved for analyses of more complete mass models. It might also be noted that 
the apparently simpler procedure of adjusting initial conditions for the forward integration 
instead of final positions for NNAM proves to be less efficient because is a smoother 
function of present positions. 



over the parameters that are adjusted — MW circular velocity, particle masses, and Cartesian 
components of present particle positions — and the parameters that are derived — redshifts, 
proper motions, and initial peculiar velocities. 

NNAM provides orbits with initial velocities that approximate the growing mode of 
departure from homogeneity, but that leaves free the magnitude of the initial velocity, which 



3.2. Physical parameters 



The measure of how well the solutions fit the constraints is the sum 




(12) 
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Ta.l)le 1: Constraints 



quantity 


catalog 


1 


2 


3 


4 


5 


Vc 


230 ± 10 km s"^ 


233 


230 


236 


236 


233 


m(MW) 


15 X X lO^mo 


12.4 


15.1 


17.8 


6.5 


12.8 


m(M31) 


30 X 2±i X lO^m© 


17.2 


24.0 


12.3 


11.5 


20.7 


m(LMC) 


0.6 X X lO^^m© 


0.92 


1.91 


0.20 


1.48 


0.93 


m(M33) 


1 X 3^1 X 10^^7710 


0.71 


0.62 


0.44 


4.6 


5.2 


m(IClO) 


0.5 X 3±^ X lO^^mo 


0.59 


0.53 


0.22 


0.20 


0.44 


m(N3109) 


1 X 3±^ X lO^^mo 


2.3 


0.54 


5.4 


0.58 


0.48 


m(N300) 


2 X 3±^ X lO^^mo 


2.2 


4.2 


2.3 


1.38 


0.79 


m(Mafrei) 


60 X 3±^ X lO^^mo 


57 


85 


87 






cz{M31) 


-300 ± 30 km s"^ 


-257 


-275 


-279 


-293 


-318 


cz{LMC) 


278 ± 10 km s^^ 


278 


276 


280 


280 


276 


cz{M33) 


-179 ± 20 km s"^ 


-196 


-192 


-219 


-227 


-206 


cz(IClO) 


-348 ± 20 km s"^ 


-349 


-345 


-339 


-320 


-319 


c;2(N3109) 


403 ± 30 km s"^ 


358 


364 


303 


392 


356 


cz(N300) 


144 ± 40 km s"^ 


136 


120 


134 


135 


135 


c^(Maffei) 


48 ± 100 km s-^ 


25 


26 


30 






D(M31) 


0.785 ± 0.1 Mpc 


0.79 


0.91 


0.67 


0.78 


0.79 


D(LMC) 


49 ± 6 kpc 


53 


55 


58 


50 


51 


D{M33) 


0.81 ±0.1 Mpc 


0.74 


0.68 


0.88 


0.75 


0.67 


D{IC10) 


0.76 ± 0.2 Mpc 


0.96 


1.21 


0.85 


1.02 


1.29 


D(N3109) 


1.3 ±0.1 Mpc 


1.37 


1.37 


1.43 


1.31 


1.37 


i:'(N300) 


2.0 ±0.4 Mpc 


2.11 


2.25 


2.06 


2.00 


2.15 


L)(Maffei) 


3.1 ± 1.0 Mpc 


3.1 


3.4 


2.7 






d±{M31) 


±10 kpc 


6.7 


4.1 


0.9 


0.5 


9.7 


di(LMC) 


±2 kpc 


0.8 


2.0 


6.4 


1.3 


5.1 


dx(M33) 


±6 kpc 


1.4 


1.3 


1.7 


1.8 


2.8 


c^^(IClO) 


±6 kpc 


0.7 


1.5 


1.5 


6.1 


0.4 


dx(N3109) 


±10 kpc 


0.3 


0.3 


2.0 


0.6 


0.4 


d±{moo) 


±0.3 Mpc 


0.06 


0.12 


0.22 


0.10 


0.05 


ci^(Maffei) 


±1 Mpc 


1.3 


0.5 


1.2 






/^a(LMC) 


2.03 ± 0.08 mas y"^ 


2.05 


2.08 


2.05 


1.99 


2.08 


/i5(LMC) 


0.44 ± 0.05 mas y"^ 


0.44 


0.44 


0.47 


0.42 


0.43 


;U„(M33) 


23 ± 6 /xas 


22 


15 


20 


22 


5 


At5(M33) 


2 ± 7 /xas y~^ 


11 


8 


4 


-4 


6 


/^a(IClO) 


— 2 ± 8 fxas 


13 


10 


9 


24 


14 


/^^(ICIO) 


20 ± 8 fias y-^ 


10 


20 


16 


-11 


-10 






25 


30 


48 


46 


61 
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corresponds to the freedom of choice of the amphtude of the growing departure from ho- 
mogeneity in an ideal fluid. Initial velocities are constrained by assigning = 50 km s~^ 
to each Cartesian components of the initial peculiar velocity of each galaxy. This choice is 
motivated by the thought that initial peculiar velocities might be expected to be comparable 
to velocities within growing protogalaxies. The solutions presented in the next section have 
reduced ~ 0.3 for the initial velocities, meaning NNAM produces smaller initial velocities 
than the adopted estimate. But this term in the total is needed because NNAM also 
produces solutions with much larger initial velocities. (The sum over ignores the fact 
that the center of mass is at rest, meaning velocities are slightly overcounted, but correcting 
this does not seem worth the trouble.) 

The cosmological parameters in the Friedmann equation ([T]) are fixed at 

n = 0.27, Ho = 70 km s"^ Mpc'^ (13) 

Limited experiments with other choices within the current uncertainties indicate the model 
results are not sensitive to these parameters. 

The second column in Table 1 lists catalog or nominal values of the parameters other 
than initial velocities in the sum. Central values of angular positions, redshifts, and 
most distances are based on NED. Many of the central values of masses are from an earlier 
study of the LG and its surroundings (Peebles, Phelps, Shaya, & TuUy 2001), but the values 
for LMC, M31 and ICIO are guesses based on luminosities. Most assigned uncertainties (jj 
in positions and redshifts are larger than measurement errors. They may be taken to be 
estimates of how far the visible galaxies might be offset from the effective centers of the 
dominant dark matter halos, or as a way to allow for the approximate nature of the model. 
The quantity d± in the table is the adopted standard deviation in each of the two orthogonal 
components of the linear distance between catalog and model angular positions (where d± is 
measured from the catalog position along the line perpendicular to the catalog direction to 
the line in the angular direction of the model). The masses are imagined to have lognormal 
distributions, that is, the difference of logarithms of model and catalog masses enters the 



square in equation (12). All other parameters are treated as normal variables. 



The first entry in the table is the circular velocity Vc of MW at our position. The 
central or nominal catalog value is taken to be intermediate between the standard and a 
possibly larger value (Reid et al. 2009), and the adopted uncertainty allows the standard at 
one standard deviation. The choice of Vc affects the model for the MW mass distribution 
(eq. [9]), but more important is the relation between heliocentric and galactocentric velocities. 
The Solar velocity relative to the local standard of rest is taken to be f/ = 11.1, V = 12.2, 
W = 7.2 km s~^ (Schonrich et al. 2010). This is fixed, not an entry in x^. 
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The central value for the mass of MW is in the range of usual estimates, with an 
allowance of twice or half this value at one nominal standard deviation. The effective center 
of the mass associated with MW is taken to be at the center of the galaxy, fixed at 8.5 kpc 
distance from the Sun. It would more consistent but overly complicated to allow an offset of 
this position along with the other galaxies. The mass of M31 is assigned a larger central value, 
as is thought to be likely, again with a factor of two allowance at one standard deviation. 
The uncertainty in the effective redshift of the mass belonging to M31 is taken to be about 
10% of the circular velocity. This seems reasonable but, as for many of the other parameter 
uncertainties, it is an intuitive guess. The distance is from the survey by McConnachie et 
al. (2005), but the adopted standard deviation is larger than the stated measurement error. 
The allowed perpendicular distance d±^ of the effective center of the mass of M31 from the 
observed center of the galaxy is 10 kpc, about the optical width of the galaxy. 

The mass of LMC is allowed a larger range because it is interesting to see whether the 
dynamics give some indication of its likely value. The distance and its uncertainty are from 
Freedman et al. (2001). van der Marel et al. (2002) describe an uncertainty of about 1 kpc 
in the position of the optical center of the LMC. The choice dj_ = 2 kpc for the offset from 
the dark matter halo thus seems justified but maybe overly optimistic. It is adopted because 
an offset much larger than this significantly shifts the LMC angular position, confusing the 
meaning of the measured proper motions. The LMC proper motions in Table 1 are from 
Kallivayalil et al. (2006), and the adopted standard deviations for are their stated errors 
(with no allowance for possible motion of the stars relative to a dark matter halo) . 

The galaxies M33 and ICIO are in this analysis because their proper motions (Brunthaler 
et al. 2005; Brunthaler et al. 2007) are important constraints. The stated uncertainties in 
the measured proper motions are treated as standard deviations. The allowed perpendicular 
offset d± is more generous than for the LMC because the angular positions are much less 
sensitive to d±. The larger allowance in the distance to ICIO seems warranted by its low 
galactic latitude (Sanna et al. 2008; Kim et al. 2009). 

The galaxy NGC3109 is included because this relatively small spiral with its scattering 
of dwarf companions is in a low density region on the edge of or just outside LG. At this 
location the redshift and distance of NGC3109 is expected to give a particularly direct 
constraint on the LG mass. The distance is from Dalcanton et al. (2009), with standard 
deviation larger than the measurement error. 

The object NGC300 is meant to represent the mass belonging to this galaxy and to 
M55. It may also represent the tidal field of the more distant galaxies in the Sculptor 
Group. To accommodate this, the allowed ranges of position and redshift are broader than 
the assignments for the nearer objects. The object Maffei is similarly meant to represent the 
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mass around IC342 and Maffei 1 and perhaps also the tidal effect of more distant mass in 
roughly the same direction. This first exploration of how the orbit of the LMC relative to 
the MW might be affected by large mass concentrations external to the LG compares models 
with and without Maffei. 



4. Models 

Table 1 lists parameters in three solutions that include the large external mass in the 
direction of the Maffei system of galaxies, and two solutions without this external mass. All 
five models fit the measured motion of LMC. The measure of fit to all the constraints 
is entered in the bottom row of Table 1. It is the sum over squares of differences of model 
and nominal values of fc, 8 masses, 7 redshifts, 3x2 proper motions, 7x3 components of 
present positions, and 7x3 initial velocities, for a total of 64 parameters that are imagined 
to have Gaussian distributions. The 30 quantities adjusted to reduce ^ire fc, the 8 masses, 
and the 7x3 components of present positions. Thus we might expect ~ 34 in a realistic 
solution, if the uncertainties o", were realistic. 

The particle orbits in comoving galactic coordinates are shown in Figures [T] and |2] They 
show two types of motion of LMC and MW. Type I, in Models 1, 2 and 5, has LMC initially 
moving up in a clockwise direction in the xz and yz projections. Type II, in Models 3 and 4, 
has LMC more directly approaching MW from above with a sharp bend at the low redshift 
end. These two general orbit types, with considerable variations in the orbits of M33 and 
ICIO, are identifiable in all solutions I have found with x^ ^ 60. The two orbit types largely 
differ by the motion of MW relative to the center of mass of the mass model — the motion 
of LMC relative to MW is quite similar, as will be discussed — but there are two possibly 
significant points to note. First, parameter adjustments in the walk to lower x^ for Type II 
solutions often end at an abrupt change of orbits and a large increase in x^. The walk in 
Type I seems to be approaching a smooth minimum that is not quite reached in acceptable 
computation time. Second, limited trials with the initial expansion factor changed from 



equation (11) to Oj = 0.2 yielded Type II solutions that are very close to what is found at 
Oi = 0.1, but did not yield Type I solutions. 

Models 1 and 2 have the smallest x^ found in this study. They have similar parameters, 
and the orbits are similar but distinctly different. Both are shown to illustrate the variations 
allowed in similar arrangements of the orbits at different apparently local minima of x^. 
Models 3 and 4 are included to show the alternative Type II behavior. Models 4 and 5 
illustrate the effect of eliminating the large external mass, at about the lowest x^ for this 
purpose. 
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Fig. 1. — Orbits in models 1 to 3 in Table 1 in a right-hand comoving coordinate system with 
X-axis along galactic I = b = and 2;-axis along b = 90° . The center of mass is at rest and the origin 
placed at the present position of the Milky Way. Labels are near initial positions at 1 -|- = 10. 




-2 2 -2 2 

X Mpc y Mpc 

Fig. 2. — Same as Figure [T] for models 4 and 5 in Table 1. 

All 5 models put Vc above 220 km s^^. This may not be significant, however, because 
the nominal catalog value is larger too, and may have biased the choice of local minima. 
Remaining to be done is a search for solutions based on a smaller catalog value of Vc- 

Since the masses are allowed considerable departures from nominal it is a positive result 
that in Models 1 and 2 the mass of M31 is larger than MW, as is usually considered likely. The 
sums of masses of M31 and MW are 3 x lO^^m© and 4 x lO^^m©, a substantial difference 
largely due to the allowed difference of model redshifts of M31. In these two models the 
components of proper motion of ICIO differ from the measurement by less than two times 
the stated error, which seems acceptable, and the other proper motions are closer to the 
measurements. The distance of ICIO in model 2 is large, but obscuration complicates this 
measurement. These models place Maffei well away from its catalog position, at c?_l = 1.3 
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and 0.5 Mpc, but this is an acceptable outcome of the idea that Maffei may represent external 
mass in the general direction of the nearest large galaxy concentration. In both models is 
consistent with what is expected from the parameter count. But since many of the nominal 
catalog standard deviations are not much better than guesses the conclusion is that the 
models are reasonable fits to what is thought to be known within the ambiguities proposed 
in Section 1312! 

Model 3 has a reasonably small but three questionable features. First, it puts LMC 
near the edge of the optical image, at d±_ = 6.4 kpc. Second, MW is more massive than M31, 
which seems unlikely though perhaps not impossible. Third, the redshift of NGC3109 is 100 
km s~^ below the catalog value, which again seems unlikely. To be checked is whether a 
more complete mass model would allow the Type II pattern of motion to better fit arguably 
reasonable constraints. 

Models 4 and 5, without the large external mass, have acceptable proper motions of 
LMC and proper motions of M33 that are arguably not unreasonable. The proper motion of 
ICIO is off by nearly four standard deviations, but this is may be acceptable for the purpose 
of modeling the motion of LMC, because at its greater distance ICIO may have been affected 
by more distant objects not in this mass model. Perhaps the greatest objection to Models 4 
and 5 is that they replace the external mass with a large mass of M33, at roughly half the 
mass of the MW, which seems unlikely. But the important issue for the present study is how 
this rearrangement of the more distant mass affects the orbit of LMC relative to MW. 

Figure |3] shows the history of angular positions of LMC and M31 for an observer fixed 
in the MW at the present position of the Sun in the galaxy. The present angular position of 
LMC is near the bottom of the figure. Model 3, with its questionable value of d±, is most 
different from the other orbits. The considerable scatter among the others includes a clear 
difference between the orbits of the two most plausible Models 1 and 2 (plotted in black). 
The orbits share some distinctive features, however. All except Model 3 passed within ~ 20° 
of the South galactic pole, in the general direction of the Magellanic stream (Mathewson & 
Ford 1984). Within the scatter these models agree with the Besla et al. (2007, Fig. 8) model 
for the motion of LMC past the pole. At redshift near unity all five orbits pass a stationary 
point in latitude near 6 = 0. At this point the longitude is near / = 90° and decreasing. This 
behavior is seen also in the more schematic mass model with cosmological initial conditions 
in P9. 

The present position of M31 is toward the left side of Figure [3] In all models M31 is 
moving to increasing galactic longitude, in the opposite direction to the motion of LMC. 
The net displacement of M31 is smallest in Models 4 and 5 (and the blue curve is so short 
it is difficult to see in the figure), as might be expected from the absence of torquing by a 
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Fig. 3. — Angular positions of LMC and M31 relative to the present solar position, for Model 1 
(solid black), 2 (dashed black), 3 (solid red), 4 (dashed red) and 5 (blue). The curves converge on 
present positions of LMC near the bottom of the figure and of M31 near the left-hand side where 
the shorter curves converge. The squares on LMC orbits mark positions at z = 0, 1, and 3. 



large external mass. The heliocentric velocity of M31 is Va = 234 km s~^, vs = —65 km 
in Model 1 and Va = 223 km s~^, vs = 8 km s~^ in Model 2. van der Marel & Guhathakurta 
(2008) argue for f q, = 78 ± 41 km s~^, vs = —38 ± 34 km s~^, which differs by nearly four 
times the stated error. Understanding this discrepancy requires further discussion of both 
approaches, and, it is to be hoped, the situation will be established by direct proper motion 
measurements in M31. 

Figure |4] gives a three-dimensional picture of the orbits of LMC relative to MW at the 
present Solar position. One sees again the considerable differences but also the similarity of 
the patterns of motion. In Model 1, LMC is 190 kpc away from MW and 300 kpc from M31 
at 1 + 2; = 10, and LMC reaches maximum distance 500 kpc from MW at z = 0.66. The 
numbers are quite similar in Model 2. Given the MW model mass, the initial and maximum 
LMC-MW separations are similar to that of a radial orbit that leaves an isolated MW at 
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high redshift and has just returned. That is, the orbit is largely dominated by MW, though 
the other galaxies are needed to account for the development of this orbit out of cosmological 
initial conditions. 



5. Conclusions 



The case that Models 1 and 2 usefully approximate the past motion of the Large Magel- 
lanic Cloud relative to the Milky Way (in the black curves in Figs, [s] and |4]) commences with 
the existence of these arguably plausible fits to the parameters in Table 1. This result was 
not guaranteed, and it is important therefore that the search has yielded one — and only one 
— fully acceptable arrangement of orbits, with the variations illustrated by the differences 
between Models 1 and 2. The values of in these two models are consistent with the counts 
of constraints and adjustable parameters within the nominal standard deviations in Table 1. 
Since most of these standard deviations are based on intuitive arguments this result should 
be read to mean that Models 1 and 2 are plausible within the arguments of reasonableness 
in Section [3j One could attempt to develop a more quantitative assessment of the meaning 
of the values of by applying the analysis to random catalogs, but that does not seem 
worthwhile; the case is better made by broader considerations. 
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One such consideration is that the nature of the orbit of LMC relative to MW is similar 
in all five models (Figs. |3]and|4]). These models all fit the measured motion of LMC, and they 
generally fit all the measured redshifts and positions in the mass model within the assigned 
uncertainties. Model 3 has problems with the larger mass of MW, the large offset of LMC 
from the optical center, and the small redshift of NGC3109, but it shares the single close 
passage of LMC and the swing to larger galactic longitude going backward in time at redshift 
about unity. Models 4 and 5 illustrate the need for three massive actors in addition to LMC 
to account for the motion of LMC out of the LMC-MW-M31 plane. (As noted in P9, the 
gravitational interaction among three galaxies, LMC, MW and M31, in an otherwise empty 
universe drives motions confined to the plane of the three galaxies.) The models without a 
large external mass solve the problem by promoting M33 to a massive actor. This is unlikely 
but a useful illustration of what the measured LMC motion requires. And it is to be noted 
that the LMC orbits relative to MW in Models 4 and 5 share the general features of the 
other solutions. The same is true of the still more schematic model in P9. 

Since the motion of LMC at low redshifts is dominated by the mass in MW, it is not 
surprising that the orbits in Figure [3] at low redshift all are similar and agree with models 
with MW alone (Besla et al. 2007) or with MW and M31 (Kallivayalil et al. 2009; Shattow 
& Loeb 2009). The swing toward larger galactic longitude at larger redshift, and the single 
close approach of LMC to MW, are less intuitive, but they have proved to be stable results 
of the measured LMC motion under the cosmological initial condition that the primeval 
peculiar velocities are growing. 

Besla et al. (2007) and Shattow & Loeb (2009) present models in which LMC has 
completed more than one orbit around MW. This is not found in the present analysis. The 
random starting orbits were designed to reach plausible solutions with two close passages, 
but the lack of discovery of examples must be balanced against the experience that the action 
method is not well adapted to capturing this case. The argument from astronomy against 
an earlier close passage is that it might be expected to have stripped away the HI (van den 
Bergh 2006). Thus Besla et al. (2010) find that the single orbit allows a model for the 
origin of the Magellanic stream by tidal interaction between the Large and Small Magellanic 
Clouds. 

The computation allows considerable freedom in the mass of LMC, and indeed the mass 
in Model 3 is one tenth that of Model 2. However, the more plausible Models 1 and 2 put the 
mass of LMC at 1 to 2 x IO^^ttiq. This is ten times the mass within 10 kpc (van der Marel 
et al. 2002), but at luminosity ~ 3 x IO^Lq a conventional dark matter halo could make 
up the difference. Indeed, Besla et al. (2010) present an independent argument for an LMC 
mass ~ 2 X lO^^m0 from the relation between luminosity and halo mass found in numerical 
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simulations of structure formation in the standard cosmology. If LMC had a massive halo 
as it approached MW then this dark matter would now be far from smoothly distributed 
as the halo merges with the MW halo. A halo merger does not seem likely to have had a 
significant effect on the motion of the stellar part of LMC as it plunges into the dominant 
MW dark halo, but a numerical analysis to check this would be interesting. 

These arguments make a good — though not definitive — case that the orbit of LMC 
relative to MW has the shape illustrated in Figures [3] and |4j An issue to be revisited with 
better mass models is that this analysis depends on the assumption that the Magellanic 
Clouds originated as clumps of matter left little disturbed by the assembly of the major 
galaxies. An alternative picture in which the Clouds are debris from a violent major merger 
at a more modest redshift (Yang & Hammer 2010) is not likely to be found by NNAM. 
This picture will have to be evaluated by other considerations, including observational and 
theoretical studies of the baryonic debris from mergers. 

The difference between the LMC orbits in the two plausible models (plotted in black) 
shows that under the cosmological initial condition there still is considerable uncertainty in 
where LMC was at redshift z = 1. This might be reduced by tighter modeling of what more 
distant galaxies were doing. 

Though models 1 and 2 are reasonable fits to the measured proper motions of M33 and 
ICIO as well as LMC, the M33 orbit does not agree with the proposal by McConnachie et 
al. (2009) that M33 passed close to M31. Since the orbits of these two galaxies are sensitive 
to the orbit of M31 a firmer case for where M33 and ICIO have been awaits a firmer case for 
the proper motion of M31. In all five models in Table 1, M31 is moving toward increasing 
longitude, and in the two most plausible cases the motion to increasing right ascension is 
faster than the van der Marel & Guhathakurta (2008) estimate. Since the transverse motion 
of M31 is expected to be more sensitive than LMC to the mass distribution outside LG 
a firmer case for the proper motion of M31 and the orbits of M33 and ICIO awaits more 
detailed mass models that more completely account for the large galaxies exterior to the 
Local Group and, equally important, include more of the nearby isolated smaller galaxies 
that largely serve as test particles. The analysis presented here shows that NNAM will be 
be well suited for analyses of these more ambitious mass models, though likely with a more 
efficient way to minimize the measure of fit and a better strategy to search for acceptable 
orbits, perhaps along the lines of Peebles et al. (2001). Most important, of course, will 
be the tighter constraints from advances in measurements of galaxy distances and proper 
motions, ground-based and from the Gaya and SIM lite satellite missions. 

This analysis is based on a picture for the mass distribution and a cosmology that have 
passed searching tests. But the cosmology and our ideas about how mass is distributed 
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around galaxies still are enormous extrapolations from what actually is well established. 
The apparently successful fit to the motions of the very nearby galaxies is a modest but 
not trivial addition to our fund of cosmological tests, and a test that can and should be 
considerably improved. 

I have benefitted from advice from Ed Shaya and Brent Tully. This research has made 
use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propul- 
sion Laboratory, California Institute of Technology, under contract with the National Aero- 
nautics and Space Administration. The numerical matrix inversion of equation (A15) is from 
Press et al. (1992). 



A. New Numerical Action Method 

This describes the NNAM solution of equation ([T]) for the adjustment of the coordinates 
of a single particle toward a solution of the equation of motion. Finding simultaneous shifts 
of all coordinates of all particles follows by the same argument. It is not presented here 
because I have found shorter computation times for iterated adjustments of coordinates of 
one particle at a time. 

Let 

2 • 2 ' 

'^n+l/2^"+l/2 ^n-l/2'^«-V2 i dtn tn+l/2 — tn-1/2 

" " TT^' ^n=——— =^n-l. — = . (Al) 

and for simplicity drop the label for the particle i whose position is being adjusted. Then 
equation ([s]) is 



(A2) 



with Ff = 0. 



The equation to be solved for the coordinate shifts 5xk,n of the particle whose orbit is 
being adjusted is 

Sk.n + / , Sk n.k'n'SXk' n' = 0, — - = Sk n;k',n' = Sk' n':k,n- (^3) 

^ ' OXk,nOXk',n' 

The nonzero derivatives are 

Sk,n;k,n+1 = ^ Sk,n;k,n-1 = ' (^4) 



Sk,n,k',n — (F^ + )5k,k' H \ rf^'^ H ;r^^k,k' ) 



Since Sk,n;k',n' = unless n' = n or else n' = n±l and k' = k, equation (A3) is 



,n;k' ,n^Xk' ,n ~\~ Sk,n;k,n~ l5Xk,n-l = 0. 

k' 

Set — )■ 71 — 1 in this equation and rearrange it to get 

e _ Sk,n-l + Sk,n-l;k> ,n-l5Xk> ,n-l + Sk,n-l;k,n-25Xk,n-2 

i^k,n—l;k,n 

This gives 5xk,n in terms of 5xk,n-i and 5a:;fc^„_2. On iterating we arrive at the form 



(A5) 



(A6) 



k" 



At n = 1 this is just 



Ak,i — 0, Bk^i-k" — Sk,k"- 



At the second time step from the start, n = 2, equation (|A7|) is 

SXk,2 = 



,l;k',lSXk\l I Sk,l-k,2 



because there is no Sxkfl- Comparing this with equation (A8) we see that 



Ak,2 — —Sk,l/ Sk,l-k,2, Bk,2;k' — — Sk,l-k' ,1 / S k,l-k,2- 



(AT) 



(AJ 



(A9) 



(AlO) 



(Aii: 



At n > 3 the result of substituting the form (|A8|) into equation (|A7|) is 

^•^k,n ~ 



Sk,n-1 + ^Sk {Ak',n-1 + ^Bk',n-l;k"SXk",l) 

k' k" 

+ Sk ,n—l;k,n—2 Bk,n-2;k"SXk",l) / Sk ,n—l;k,n- 



(A12) 



k" 



This shows that at 3 < n < n^^ the constants in equation (|A8|) are 

A-k,n = 
Bk,n;k" 



Sk,n-1 + ^fc/ Sk^n-l;k' ,n-lA.k' ,n-l + 'S'fc,n-l;fc,n-2^fc,n-2 

^k,n~l;k,n 

Xlfc' Sk,n~l;k' ,n-lBk' ,n-l;k" + Sk,n-l;k,n-2Bk,n~2;k" 



Sk,n—l;k,n 



(A13) 



Equations (All) and (A13) give the A^n and B^^n^k''^ 2 < n < rix, in terms of the 



given derivatives of the action. Then equation (AS) is Sn^ — 3 equations for the Sxk^n, at 
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2 < n < n^, in terms of the 6xk,i- We get three more, which fix the Sxk^i, by setting n = 
in equation (A6) and recalhng that 5xkn^+i = 0: 



k' 

= Sk^n^ + Sk,n^-k' ,11:,: ^fc',nx + By n:,;k"^Xk" 



rix-l 



(AM) 



k" 



^fc,n^-l + Bk n:,:-l-k"^Xk" ,1 



or 



(A15) 



where 



-^A; Sk,nx ~^ Sk,n:,:;k' ,na:-^k' ,nx ~^ Sk,nx\k,nx — l-^k,n:,: — l: 

k' 

Tk,k" = Sk,nx;k',nxBk',nx;k" + Sk,nx;k,nx-lBk,nx-l;k" ■ (Al6) 

fc' 

This is three equations for the three Sxk,i- The rest of the Sxk,n then follow from equations 



(A8) and (A13). 



For Np particles and rix time steps the computation time to get the Sxi^k,n for all particles 
scales as N^N^. In previous applications of the numerical action method (P9 and references 



therein) equation (A3) is solved by matrix inversion, an operation that scales as N!:n^ for 



relaxation of one particle orbit at a time. The NNAM operation scales in the same way 
as a conventional numerical integration, but the numerical prefactor is considerably larger 
because it takes many iterations of the coordinate shifts Sxi^k,n to drive the Si^k,n to zero. 
NNAM certainly will not replace conventional forward numerical integration for simulations 
of the growth of cosmic structure. But the application presented here shows why NNAM 
will be useful for analyses of the dynamics of the nearby galaxies. 
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